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Abstract 



This work approaches human chromosome mapping by developing 
algorithms for ordering markers associated with radiation hybrid data. 
Motivated by recent work of Boehnke et al. [1], we formulate the ordering 
problem by developing stochastic spin models to search for minimum-break 
marker configurations. As a particular application, the methods developed 
are applied to 14 human chromosome-21 markers tested by Cox et al. [2]. 
The methods generate configurations consistent with the best found by 
others. Additionally, we find that the set of low-lying configurations is 
described by a Markov-like ordering probability distribution. The 
distribution displays cluster correlations reflecting closely linked loci. 
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I. INTRODUCTION 



The use of data from radiation hybrid (RH) experiments has become a 
useful method for fine structure mapping of human chromosomes. Based on 
methods described by Goss and Harris [3, 4], Cox et al. [2] and Burmeister 
et al. [5] have developed the technique in detail so that results from their 
experiments provide material for ordering DNA markers on human chromo- 
somes. 

The basic strategy employed in radiation hybrid mapping (fully described 
in Cox et al. [2]) entails irradiating a rodent-human somatic cell hybrid, that 
contains a particular human chromosome, with a lethal dose of x-rays. This 
will cause the chromosomes to break into several fragments. After fusion with 
HPRT-deficient rodent cell lines, only the fused cells, containing both the X- 
irradiated cells and the normal rodent cells, will survive if grown in HAT 
medium. Each hybrid clone arising from this fusion will contain a unique 
set of fragments from the original human chromosome. Each clone can then 
be typed for a set of human DNA markers (equivalently, loci) known to be 
on that human chromosome. Based on the assumption that tightly linked 
markers are unlikely to be broken apart by the radiation, markers close to 
one another will show a correlated pattern of retention in the clones; whereas, 
distant markers will be retained in a relatively independent manner. 

Several methods for ordering markers have been developed using results 
from RH experiments [including both parametric (Cox et al. [2]; Boehnke et 
al. [1]) and non-parametric methods ( Boehnke et al. [1]; Falk [6]; Bishop 
and Crockford [7]; Weeks et al. [8]). In particular, Boehnke et al. [1] used 
a mathematical quantity associated with the number of breaks and then 
used optimization techniques to minimize that quantity. One optimization 
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technique involved a simulated annealing search for configurations associated 
with minimal numbers of breaks. 

We set out to study and understand the work of Boehnke et al., and we 
developed a formulation in the context of stochastic spin systems (Falk [9]). 
It may be useful to point out some differences in the formulations. 

Boehnke et al. use block inversions of a given marker order and compare 
the old and new orders with respect to "obligate" breaks. They then apply 
simulated annealing techniques and decide, at each time step, whether to 
retain the original order or transition to the new. If a transition would result 
in a smaller number of breaks, the transition is made with probability one. 
If a transition would not decrease the number of breaks, then the transition 
probability is less than one, and that transition probability systematically 
decreases over time. 

In our study we implement three algorithms which incorporate the num- 
ber of "breaks". The three algorithms are three stochastic spin models. 
These, too, are designed to search for configurations with small numbers of 
breaks. A probability is constructed to determine whether or not to retain 
the current order or transition to the new. The probabilities are set up so as 
to bias the decision towards transitions to configurations with fewer breaks; 
however, at a given step, a possible transition leading to a smaller number 
of breaks, will not necessarily be realized. 

The spin language provides mathematically intuitive expressions for the 
number of breaks. Those expressions contain products of adjacent spin vari- 
ables, and calculations involving breaks are easily presented in spin notation. 
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For those seeking a rigorous mathematical setting, we remark that Liggett 
[10] has treated stochastic spin models and related models from biology, 
physics and econonomics. Liggett 's book contains an extensive bibliography 
and guides the reader to survey papers by Griffeath, Durrett, Stroock, Hol- 
ley, and others. 

II. METHODS 

We are considering M clones, each of which is tested for the presence (or 
absence) of N different DNA markers. It is convenient to represent the clones 
as M rows, each with N sites. Thus, one pictures a two-dimensional array of 
M rows, N columns. 

1 2 N 

1 

2 



M 

Assign a particular DNA marker, labelled fj, (j = 1, 2, TV), to each 
column. Associate a variable ("a spin") Sij with the site at row i, column j. 
If the marker fj is present at site (i, j), take s^- = +1; if fj is not present 
at site take = —1. If one is uncertain as to whether a marker is 

present at site (i, j), we take s^- to be unknown, and we deal with such sites 
in a manner to be specified subsequently. 

For (j = 1,2, N) the marker fj was arbitrarily assigned to column j. 
But any of the N\ assignments would be possible. A criterion is needed for 
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judging the "goodness" of an assignment. 

Considering that the DNA markers being tested for lie on a particular 
human chromosome, those markers which are tightly linked are likely to ap- 
pear together or not appear in each clone. Therefore, it is reasonable to 
seek those marker assignments which reflect such correlation. For that pur- 
pose we say that a "break" exists between sites and + 1) whenever 
SijSij + i = —1. The strategy is to minimize the total number of breaks. Note 
that in the two-dimensional array of spins, a break refers only to horizontal, 
nearest-neighbor spin pairs. 

Here we devise and test several algorithms which monitor the total num- 
ber of breaks while selectively permuting column labels. The algorithms 
attempt to explore the vast configuration space in the manner of a "random 
walk" , biased toward configurations having a reduced number of breaks. No- 
tice, for N = 14 there are N\ = 87178291200 permutations of the numbers 
(1, 2, 14). Thus, in the spirit of the travelling salesman problem, simulated 
annealing techniques are used (Kirkpatrick et al [11], Press, et al [12]). 

Nearest-neighbor transposition algorithm 

(0) Start with a specified configuration {s} of the variables s^. 

(1) Select a number j at random from the set {1, 2, iV}. 

(2) If j 7^ N, compute the total number of breaks between columns j — 1, j 
and between columns j + 1, j + 2. Denote that number by Bj(l): 
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^•(i) = Eta - M T qj + C 1 - ^-i) % +l 9J+2 ] (2.1) 



=1 



for j G {1, 2, TV — 1}. The Kronecker delta contained in (1 — Sij) is inserted 
to handle "end effects" since column 1 has no left neighbor and N has no 
right neighbor. 

(3) Repeat (2) with the spin values and Sjj+i interchanged for i = 
1, 2, M, and denote the resulting sum by Bj(2) instead of Bj(l). 

(4) Then compute Bj, the change in the number of breaks resulting from 
the interchange of the columns of spin values s^- and Sjj+i. 



Bj = B,(2) - B 3 {1) for j £ {1, 2, N — 1}. (2.2) 

(5) Interchange DNA marker column assignments and the associated spin 
values for columns j and j + 1 with probability. 



= eXP( -^ /M) , for, e{ 1.2 V. ,)• 



for j = JV. (2.3) 



Do not interchange DNA marker column assignments and the associated spin 
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values for columns j and j + 1 with probability 



W 2 j = 1 - Wly 



(2.4) 



Note: 



w 2j = 



exp (+/3Bj/M) 
2 cosh (PBj/M) 



for je {1,2, ...,7V- 1} 



= 1 



for j = iV. 



(2.5) 



The "inverse temperature" parameter (3 ( < (3 < oo ) is allowed 
to increase in an empirically determined manner as the algorithm is imple- 
mented. Clearly very large /3 strongly favors transitions reducing the number 
of breaks; whereas, a small, positive value of (3 makes transitions to increase 
or to decrease the number of breaks almost equally likely. Why not just 
take (3 large at the outset? The answer is (Kirkpatrick et al [11], Press et al 
[12]) that setting (3 large early in the calculation may cause the algorithm to 
get "trapped" in a local minimum before performing a significant search of 
configuration space. 

Thus, the initially chosen value for (3 and the manner of increasing that 
value constitute the "art" of simulated annealing. 

(6) Return to (1) and repeat the procedure starting with the current con- 
figuration {s'} of the variables s'^. The procedure may be terminated after 
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a specified number of iterations and/or when the total number of breaks 



M N-l i _ „ „ 

g = EE f qJ+1 (2-6) 

5=1 J=l Z 

has realized acceptably small values. For each run one retains the configura- 
tions associated with the smallest values of B. 

If /3 were fixed, then the above procedure would define a finite-state, 
discrete-time Markov chain with transition probability p({s'} \ {s}) from a 
configuration {s} to configuration {«'}. Explicitly 



AT 



p(WIW) = ^E 

JV 7=1 



'T + S^S ?J+ l\ / 1 + So 



n 

9=1 



AT 

x n 

k=l;k^j,j+l 



S qk S qk 



+ 



M N 

nn 

g=l fe=l 



1 + s' t s 



qk b qk 



w 2j 



(2.7) 



Two-column transposition algorithm 



A natural extension of the nearest-neighbor transposition algorithm in- 
volves columns k,j with k > j + 1. Remove two numbers at random from 
the set {1,2,..., N}. Denote the smaller number by j and the larger by k. If 
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k — j + 1, follow the previously described nearest-neighbor algorithm, start- 
ing with step (2). If k > j + 1, proceed as follows: 

Consider the total number of breaks between columns j — l,j;j,j + 1; 
k — 1, k; k, k + 1. Denote that number by Bj k (l), where 



M r 



%(i) = £ (i-^-o 



9=1 



^ o H 1 ~ d W r 



(2i 



After interchanging spin values and for i = 1, 2, M, the number 
of breaks is 



M 



9=1 



;i - ^,i) 



^ '^ij — l^qk _|_ 1 SqkSq,j-\-l 



+ 



+ (1 



(2.9) 



Then the change (in the number of breaks) resulting from the interchange 
is denoted by Bj k , where 



B jk = B jk (2)-B jk (l). (2.10) 
9 



Now interchange DNA marker column assignments and the associated 
spin values for columns j, k for k > j + 1 with probability 



Wijk = 



exp (-/3B jk /M) 
2 cosh {(3B jk /M) 



for je {1,2, ...,N -2}. 



(2.11) 



Do not interchange DNA marker column assignments and the associated 
spin values for columns j and k with probability 



Block-flip algorithm 

As in the preceding algorithm, remove two numbers at random from the 
set {1,2, ...,N}. Denote the smaller number by j and the larger by k. If 
k — j + 1, follow the previously described nearest-neighbor algorithm, start- 
ing with step (2). If k > j + 1, proceed as follows: 



W 2 jk = 1 - Wljk- 



(2.12) 
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The block consists of columns j,...,k. Before flipping the block, the num- 
ber of breaks involving columns j — 1, j and columns k, k + 1 is 



B Uoc k{l) 



M 

E 



1 s q,j-l s qj 



+ (1 - 5 k>N ) 



1 — Sq,kSq,k+l 



(2.13) 



After flipping the block, the number of breaks is 



M 



B u ock{2) = ^ {1 _ 5 . i) 

q=l 1 



1 Sq,j — lSqk 



+ a-s. 



1 SnjS, 



k,N) 



qj b q,k + l 



(2.14) 



Then the change in the number of breaks is denoted by B^ ock ^ where 



■q block -q block /q\ -q 62ocfc/-p 

jk jk V / j k V J 



(2.15) 



Now flip the block with probability 
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w 



block 



exp {-(3B$ ock /M) 
2 cosh ((3B$ ock /M) 



for je {1,2, ...,7V -2}. (2.16) 



Do not flip the block with probability 



w. 



block 
2jk 



= 1 — W 



ljfe 



(2.17) 



The preceding method is similar to the block inversion algorithm used by 
Boehnke et al., but our transition probabilities differ from theirs. 

Unknown site content 

In any clone one may have strings of one or more sites where the DNA 
marker associated with each site in a string is unknown. Thus, the spin val- 
ues are unknown for the string. In the above algorithms such unknowns are 
dealt with in the following way. 

Consider the case where the left and right ends of a string connect re- 
spectively to known spin sj e / t and to known spin s rig ht- If Sz e /t = Sright, then 
all spins in the string are replaced by s rig ht- If s/ e /t 7^ Sright, then a fair coin 
toss is simulated. If the coin shows head (tail), then all spins in the string 
are replace by +1 (—1). 
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If a string (of unknowns) contains an end spin, then replace all spins in 
the string by the value of the connecting spin (si e f t or s right), as appropriate. 

With all spin values now specified, the number of breaks can be calculated 
for any of the above algorithms, and the relevant transition probability can 
be evaluated. The simulated annealing proceeds one step. After that step, 
all of the previously unknown spin values are again regarded as unknown. 
[Note: Due to a possible column interchange or block flip, those unknown 
spins, which previously belonged to particular strings, may now belong to 
different strings.] The above prescription for replacing unknown spin values 
by +1 or —1 is now repeated, and the transition probabilities are reevaluated. 
The simulated annealing proceeds another step, etc. 

III. APPLICATION 

As an example, consider the data presented by Cox et al. [2] relating 
to 14 markers on chromosome 21 that were tested in 99 RH clones. The 14 
markers are given in Table 1. These are the same data used by Boehnke et 
al. [1] and presented in detail in their Table 1. For our algorithms an entry 
of 1 in their table becomes +1, becomes -1 and a "?" remains unknown 
and is treated at each step as described above. 

All three algorithms were applied to the data in the 99 x 14 matrix for 
200,000 iterations. Initial values of (3 and incremental steps for increasing 
f3 were chosen. This produced a set of permutations with acceptably small 
values of B, the total number of breaks for a particular configuration. For 
each run, a ranked set of marker permutations with the smallest values of B 
was retained. 
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Table 2 lists the 2 distinct "best" orders found by each algorithm in a 
representative run. As can be seen, the first algorithm, where two nearest- 
neighbor columns are transposed, does not result in permutations with val- 
ues of B that are as low as those reached by algorithms 2 and 3. Although, 
in principle, the nearest-neighbor transposition should allow for visiting all 
possible permutations of the columns, in practice, such exploration is not ef- 
ficiently accomplished here. The large configuration space, and the empirical 
nature of selecting f3 to implement simulated annealing provide a computa- 
tional challenge for the first algorithm. The other two algorithms display 
improved ability to achieve low B values with the chosen set of parameters. 
Additionally these algorithms reach the same optimal order as that attained 
by Boehnke et al. (see their Table 2), with the same number of breaks. The 
only difference is that we have retained all 14 markers; whereas they com- 
bined markers S12 and Sill , since the latter markers were indistinguishable 
in the data matrix. Hence for each marker order in their Table 2 we would 
have 2 orders. 
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Although algorithms such as these do not assure that the marker order 
with the smallest number of breaks is the correct order, inspection of a set of 
low-lying states leads to some useful information about the stability of clus- 
ters of markers that retain their local ordering. For example, consider the 
set of unique permutations representing the 24 "best" orders obtained from 
a series of runs of the three algorithms (Table 3). We see, e.g., that S47 and 
SOD are nearest neighbors in all 24 permutations and appear as the last two 
markers in 16 permutations. Similarly, the triplet S46-S4-S52 is preserved in 
21 permutations and of these, positions 2,3 and 4 contain S46, S4 and S52 
respectively 15 times. Based on observations such as these, we looked for a 
general ordering property associated with sets of low-lying configurations. 

Markovian-like property of low-lying configurations 

Let the DNA marker at site j be denoted by /_,-, where fj is a member 
of the set {S16, S48, S46, S4, S52,...}. A configuration of sites is denoted by 
the ordered N-tuple (/i, f 2 , /at). We have used underlining to distinguish 
a DNA marker such as Sll from a spin variable such as Sn. 

For any configuration of sites, one can define strings (fi,f i+ i,... fi +m ) 
with % = 1, N; < m < N - i. 

Consider a collection C of distinct configurations. In the collection the 
probability of a string, (f h f i+m ) is denoted by P i: ..., i+m (f h f i+ i, f i+m ). 

Now consider the Markovian-like approximation 
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Pi,...,i+m(fi, fi+1, fi+m) ~ Pi,i+l,i+2(fi \ fi+1, fi+2)Pi+\,i+2,i+z{fi+l \ fi+2ifi+^)--- 
■■■Pi+m—2,i+m—l,i+m(.fi+m—2 \ fi+m—li fi+m) Pi+m—l,i+m(fi+m—li fi+m) (*) 

fori = l,...,iV-2 ; 2 < m < N - i. 

Example 1. 

Consider the collection of 24 distinct low-lying configurations given in 
Table 3. Let (the number of markers) iV = 14. Look at the cluster (f e = 
F, f 7 = GJ 8 = H, f 9 = /, /io = J) where F denotes the DNA marker SU, G 
denotes SI, H denotes S18, / denotes S8, and J denotes APP. This marker 
assignment corresponds to the ordering of the first row of Table 3. 

For the above collection we find the frequency 



-^6, 7,8,9,10 (-^ ) G, H, I, J) = 7/24 (3.1) 
and we also find the frequencies 
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P 6t7t8 (F\G,H) = 11/11 (3.2) 

P 7fit9 (G\H,I) = 7/11 (3.3) 

P 8t9 , 10 (H\I,J) = 11/11 (3.4) 

P 9 , 10 (/,J) = 11/24, (3.5) 



so the approximation (*) with i = 6,m = 4 is satisfied by the observed fre- 
quencies. 

Similarly, for the above configurations P7 t 8,9,io(G, I, J) = 7/24 and 

Pw, 1Q {G, H, I, J) » P 7 , 8 , 9 , 10 (G | if, I)P 8 , 9 , W ( H I J > W,io(f , J) 
= (7/ll)(ll/ll)(ll/24) 

= 7/24. (3.6) 



Example 2. 

Consider the same collection of 24 distinct low-lying configurations used 
in example 1. Look at the cluster (/ 2 = B, f 3 — C, / 4 = D, f 5 — E) where 
B denotes the DNA marker S48, C denotes S46, D denotes S4, E denotes S52. 
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For the above collection we find the frequency 



P 2 ,s A5 (B,C,D,E) = 12/24, (3.7) 



and we also find the frequencies 



P 2M (B\C,D) = 12/15 (3.8) 
P 3A5 (C\D,E) = 15/16 (3.9) 
P 4 , 5 (D,E) = 16/24, (3.10) 



again, the approximation (*) is satisfied by the observed frequencies. 
However, since 

P 2 AB\C) = 13/17 (3.11) 

Pz A {C\D) = 15/18 (3.12) 

P 4 , 5 (D | E) = 16/16 (3.13) 

P 2 , 3A5 (B,C,D,E) = 12/24, (3.14) 

the frequencies do not satisfy the standard Markovian approximation 

18 



P 2 , 3A5 (B,C,D,E) ^ P 2i3 (B | C)P 3A (C | D)P 4 , 5 (D \ E)P 5 (E) 
12/24 ^ (13/17)(15/18)(16/16)(16/24) 
0.50 ^ 0.425 (3.15) 

A referee has kindly pointed out that in the context of Markov random 
fields (Spitzer [13]) one could perhaps find a rigorous basis for the observed 
Markov-like property. We thank the referee for this stimulating suggestion. 

IV. DISCUSSION 

The implementation of experimental techniques using RH data provides a 
bridge between chromosome mapping data generated from families and phys- 
ical mapping data. RH experiments allow for the relative ordering of genetic 
markers that are too closely spaced to be detected by family linkage analy- 
sis. Additionally, it is not necessary to have polymorphic markers in order 
to generate useful information. Although not providing the level of resolu- 
tion of physical mapping, RH mapping can complement and confirm results 
generated by pulse field gel electrophoresis. 
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Boehnke et al. [1] have presented a full discussion of the advantages and 
disadvantages of parametric and non-parametric ordering algorithms. As 
they point out, algorithms that search for minimum break configurations do 
not allow for estimates of distances between markers, nor do they provide rel- 
ative likelihoods for one marker order over another. However, as the present 
work and the work of Boehnke et al. demonstrate, retention and inspection 
of sets of low-lying configurations yield important insight relating to the cor- 
relations of markers. 

In our study we were interested in learning what properties might be 
present in a set of configurations with relatively few breaks. It became obvi- 
ous that the set of low-lying configurations showed the clustering of particular 
markers. That was reassuring, since in complex optimization problems of the 
travelling salesman variety, one typically ends up with a set of low-lying con- 
figurations and never discovers the true absolute minimum. 

One striking feature which we discovered here is that the set of low-lying 
configurations is described by a Markov-like probability distribution. That 
distribution contains all of the observed clustering of the markers. If one en- 
larges the set of low-lying configurations to include configurations with larger 
and larger numbers of breaks, the validity of the Markov-like approximation 
deteriorates. We are not in a position to judge whether the relation (*) is 
necessarily a deep or broad result, but we regard it as interesting. 
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Table 1. 



Fourteen Chromosome-21 Markers Used in the Examples. 



APP SI S4 S8 Sll S12 S16 S18 S46 S47 S47 S52 Sill SOD 



Note: All numbered loci have a prefix of D21; APP = amyloid precursor; 
SOD = superoxide dismutase. 
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Table 2. 



Minimum-Break Orders for Representative Simulations. 

Algorithm Marker Order 

| Breaks 

1 205 APP S8 SI Sll S16 S4 S12 S18 S48 S46 S52 Sill S47 SOD 

1 207 APP S8 SI Sll S16 S4 S12 S18 S46 S48 S52 Sill S47 SOD 

2 123 S16 S48 S46 S4 S52 Sll SI S18 S8 APP Sill S12 S47 SOD 

2 123 S16 S48 S46 S4 S52 Sll SI S18 S8 APP S12 Sill S47 SOD 

3 123 S16 S48 S46 S4 S52 Sll SI S18 S8 APP Sill S12 S47 SOD 
3 123 S16 S48 S46 S4 S52 Sll SI S18 S8 APP S12 Sill S47 SOD 

See note for Table 1. 

Table 3. 

24 Unique Marker Orders With Relatively Small Numbers of Breaks 
Breaks Marker Order 
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123 


S16 


S48 


S46 


S4 


S52 


Sll 


SI 


S18 


S8 


APP 


Sill 


S12 


S47 


SOD 


125 


S16 


S48 


S46 


S4 


S52 


Sll 


SI 


S18 


APP 


S8 


S12 


Sill 


S47 


SOD 


126 


S48 


S16 


S46 


S4 


S52 


Sll 


SI 


S18 


S8 


APP 


S12 


Sill 


S47 


SOD 


126 


S48 


S16 


S46 


S4 


S52 


Sll 


SI 


S18 


S8 


APP 


Sill 


S12 


S47 


SOD 


127 


Sll 


SI 


S16 


S48 


S46 


S4 


S52 


S18 


S8 


APP 


Sill 


S12 


S47 


SOD 


127 


S16 


S48 


S46 


S4 


S52 


Sll 


SI 


S18 


S8 


APP 


S12 


Sill 


S47 


SOD 


127 


S46 


S48 


S16 


S4 


S52 


Sll 


SI 


S18 


S8 


APP 


S12 


Sill 


S47 


SOD 


128 


Sll 


SI 


S52 


S4 


S46 


S48 


S16 


S18 


S8 


APP 


S12 


Sill 


S47 


SOD 


128 


Sll 


SI 


S52 


S4 


S46 


S48 


S16 


S18 


S8 


APP 


Sill 


S12 


S47 


SOD 


128 


S16 


S48 


S46 


S4 


S52 


SOD 


S47 


S12 


Sill 


APP 


S8 


S18 


SI 


Sll 


128 


S16 


S48 


S46 


S4 


S52 


SOD 


S47 


Sill 


S12 


APP 


S8 


S18 


SI 


Sll 


129 


Sll 


SI 


S16 


S48 


S46 


S4 


S52 


S18 


S8 


APP 


S12 


Sill 
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See note for Table 1. 
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